/*
代码已经经过了优化，请不要使用任何容器替换参数，务必使用线性内存，使用指针
在不开优化的情况下，使用容器与此方案差距不大，但在O2及以上优化级别差距显著
*/

#include "sampleentropy.h"
#include <math.h>

namespace entropys {

// Used by SampleEntropy()
static inline double step(double *X, int N, double r, int m)
{
    int B = 0;
    for (int i = 0; i <= N - m; i++) {

        for (int j = 0; j <= N - m; j++) {
            if (i != j) {
                double D = fabs(X[i] - X[j]);
                for (int k = 1; k < m; k++) {
                    double t = fabs(X[i + k] - X[j + k]);
                    if (D < t)
                        D = t;
                }
                if (D <= r)
                    B++;
            }
        }
    }
    return 1.0 * B / (N - m) / (N - m + 1);
}

double SampleEntropy(double *X, int N, double r, int m)
{
    double s1 = step(X, N, r, m);
    if (s1 == 0)
        return 0;
    double s2 = step(X, N, r, m + 1);
    if (s2 == 0)
        return 0;
    return -log(s2 / s1);
}

double FastSampleEntropy(double *X, int N, double r, int m)
{
    int B1 = 0, B2 = 0;
    int LoopsSub1 = N - m;
    for (int i = 0; i <= LoopsSub1; i++) {
        for (int j = 0; j <= LoopsSub1; j++) {
            if (i != j) {
                double D = fabs(X[i] - X[j]);
                for (int k = 1; k < m; k++) {
                    double t = fabs(X[i + k] - X[j + k]);
                    if (D < t)
                        D = t;
                }
                if (D <= r)
                    B1++;
                if (i != LoopsSub1 && j != LoopsSub1) {
                    double t = fabs(X[i + m] - X[j + m]);
                    if (D < t)
                        D = t;
                    if (D <= r)
                        B2++;
                }
            }
        }
    }
    double s1 = 1.0 * B1 / (N - m) / (N - m + 1);
    double s2 = 1.0 * B2 / (N - m - 1) / (N - m);
    if (s1 == 0 || s2 == 0)
        return 0;
    return -log(s2 / s1);
}

double FastSampleEntropy_ai(double *X, int N, double r, int m)
{
    int B1 = 0, B2 = 0;
    int LoopsSub1 = N - m + 1;
    for (int i = 0; i < LoopsSub1; i++) {
        for (int j = i + 1; j < LoopsSub1; j++) {
            double D = fabs(X[i] - X[j]);
            for (int k = 1; k < m; k++) {
                double diff = fabs(X[i + k] - X[j + k]);
                if (diff > D) {
                    D = diff;
                }
            }
            if (D <= r) {
                B1++;
            }
            if (j < N - m) {
                double diff = fabs(X[i + m] - X[j + m]);
                if (diff > D) {
                    D = diff;
                }
                if (D <= r) {
                    B2++;
                }
            }
        }
    }
    double s1 = 1.0 * B1 / ((N - m + 1) * (N - m) / 2);
    double s2 = 1.0 * B2 / ((N - m) * (N - m - 1) / 2);
    if (s1 == 0 || s2 == 0) {
        return 0;
    }
    return -log(s2 / s1);
}

double FastSampleEntropy_m2(double *X, int N, double r)
{
    int B1 = 0, B2 = 0;
    int LoopsSub1 = N - 2;
    for (int i = 0; i <= LoopsSub1; i++) {
        for (int j = 0; j <= LoopsSub1; j++) {
            if (i != j) {
                double D = fabs(X[i] - X[j]);
                double t = fabs(X[i + 1] - X[j + 1]);
                if (D < t)
                    D = t;
                if (D <= r)
                    B1++;
                if (i != LoopsSub1 && j != LoopsSub1) {
                    double t = fabs(X[i + 2] - X[j + 2]);
                    if (D < t)
                        D = t;
                    if (D <= r)
                        B2++;
                }
            }
        }
    }
    double s1 = 1.0 * B1 / (N - 2) / (N - 1);
    double s2 = 1.0 * B2 / (N - 3) / (N - 2);
    if (s1 == 0 || s2 == 0)
        return 0;
    return -log(s2 / s1);
}

}
